pacman::p_load(tidyverse, fishualize, patchwork, ggcorrplot)
path1      <-
    '~/Git/henriquelaureano.github.io/THESIS/code/coefs/coefs1to100/'
path2      <-
    '~/Git/henriquelaureano.github.io/THESIS/code/coefs/coefs101to200/'
path3      <-
    '~/Git/henriquelaureano.github.io/THESIS/code/coefs/coefs201to300/'
fixednames <- c('beta1', 'beta2', 'gama1', 'gama2', 'w1', 'w2')
varnames   <- c('logs2_1', 'logs2_2', 'logs2_3', 'logs2_4')
cornames   <- c(
    'rhoZ12', 'rhoZ13', 'rhoZ14', 'rhoZ23', 'rhoZ24', 'rhoZ34'
)
metnames   <- c('conv', 'mll')

vdata <- function(data, path)
{
    out <-
        tibble::as_tibble(read.table(paste0(path, data)))%>%
        dplyr::select(-V1)
    type <- substring(data, 6, 7)    
    switch(type,
           v1={ out <-
                    out%>%'colnames<-'(c(
                              fixednames,
                              varnames[1:2], cornames[1], metnames)) },
           v2={ out <-
                    out%>%'colnames<-'(c(
                              fixednames,
                              varnames[3:4], cornames[6], metnames)) },
           v3={ out <-
                    out%>%'colnames<-'(c(
                              fixednames,
                              varnames, cornames[c(1, 6)], metnames)) },
           v4={ out <-
                    out%>%'colnames<-'(c(
                              fixednames,
                              varnames, cornames, metnames)) },
           stop('Invalid type'))
    return(out)
}

v1data <- v2data <- v3data <- v4data <- vector('list', 18)
label  <- numeric(18)
l      <- 1
for (i in 1:3)
{
    for (j in 1:2)
    {
        for (k in 1:3)
        {
            label[l]    <- paste0('cs', i, '_', j, '_', k)
            v1data[[l]] <-
                dplyr::bind_rows(
                           vdata(paste0('coefsv1_', label[l], '.txt'),
                                 path1),
                           vdata(paste0('coefsv1_', label[l], '.txt'),
                                 path2),
                           vdata(paste0('coefsv1_', label[l], '.txt'),
                                 path3)
                       )
            v2data[[l]] <-
                dplyr::bind_rows(
                           vdata(paste0('coefsv2_', label[l], '.txt'),
                                 path1),
                           vdata(paste0('coefsv2_', label[l], '.txt'),
                                 path2)
                       )
            v3data[[l]] <-
                dplyr::bind_rows(
                           vdata(paste0('coefsv3_', label[l], '.txt'),
                                 path1),
                           vdata(paste0('coefsv3_', label[l], '.txt'),
                                 path2)
                       )
            v4data[[l]] <-
                dplyr::bind_rows(
                           vdata(paste0('coefsv4_', label[l], '.txt'),
                                 path1),
                           vdata(paste0('coefsv4_', label[l], '.txt'),
                                 path2)
                       )
            l <- l+1
        }
    }
}
## str(v1data);length(v1data)
## str(v2data);length(v2data)
## str(v3data);length(v3data)
## str(v4data);length(v4data)

cif1    <- setNames(c( 3,  2.6, 2.5, 4, 5, 10), fixednames)
cif2    <- setNames(c(-2, -1.5, 1, 1.5, 3, 4), fixednames)
logvars <- setNames(log(c(1, 0.6, 0.7, 0.9)), varnames)
rhoZs   <- setNames(atanh(c(0.1, -0.5, 0.3, 0.3, -0.4, 0.2)), cornames)

vtrue <- function(data.cif1, data.cif2)
{
    out <- rep(c(replicate(data.cif1, n=3, simplify=FALSE),
                 replicate(data.cif2, n=3, simplify=FALSE)), 3)
    return(out)
}
v1cif1 <- c(cif1, logvars[1:2], rhoZs[1])
v1cif2 <- c(cif2, logvars[1:2], rhoZs[1])
v1true <- vtrue(v1cif1, v1cif2)
v2cif1 <- c(cif1, logvars[3:4], rhoZs[6])
v2cif2 <- c(cif2, logvars[3:4], rhoZs[6])
v2true <- vtrue(v2cif1, v2cif2)
v3cif1 <- c(cif1, logvars, rhoZs[c(1, 6)])
v3cif2 <- c(cif2, logvars, rhoZs[c(1, 6)])
v3true <- vtrue(v3cif1, v3cif2)
v4cif1 <- c(cif1, logvars, rhoZs)
v4cif2 <- c(cif2, logvars, rhoZs)
v4true <- vtrue(v4cif1, v4cif2)

BIAS QUANTILE

cerror <- function(data, true, data_label) {
    coefs <- data%>%dplyr::select(-c(conv, mll))
    error <- coefs
    for (i in seq(ncol(coefs))) { error[i] <- true[i]-coefs[i] }
    out <- error%>%
        dplyr::summarize_all(mean)%>%
        dplyr::add_row(
                   error%>%
                   dplyr::summarize_all(quantile, c(0.025, 0.975)))%>%
        dplyr::add_row(
                   error%>%
                   dplyr::summarize_all(sd))%>%
        dplyr::mutate(label=c('mean', 'q025', 'q975', 'sd'),
                      conv=1-mean(data$conv),
                      nr=nrow(data), 
                      data=data_label)
    return(out)
}

label <- sub('cs1_', 'cs02-',  label)
label <- sub('cs2_', 'cs05-',  label)
label <- sub('cs3_', 'cs10-', label)

label <- sub('1_', 'high-', label)
label <- sub('2_', '-low-',  label)

label <- sub('-1$', '-05k',  label)
label <- sub('-2$', '-30k', label)
label <- sub('-3$', '-60k', label)

v1error <- v2error <- v3error <- v4error <- NULL
for (i in seq(18))
{
    v1error <-
        v1error%>%
        dplyr::bind_rows(cerror(v1data[[i]], v1true[[i]], label[i]))
    v2error <-
        v2error%>%
        dplyr::bind_rows(cerror(v2data[[i]], v2true[[i]], label[i]))
    v3error <-
        v3error%>%
        dplyr::bind_rows(cerror(v3data[[i]], v3true[[i]], label[i]))
    v4error <-
        v4error%>%
        dplyr::bind_rows(cerror(v4data[[i]], v4true[[i]], label[i]))
}

## v1error
## v2error
## v3error
## v4error

inside <- function(par, errordata)
{
    out <-
        errordata%>%
        dplyr::select(par, label, conv, nr, data)%>%
        tidyr::pivot_wider(names_from=label, values_from=par)
    out$data <-
        forcats::fct_relevel(forcats::as_factor(out$data), rev)
    return(out)
}
biasformat <- function(par, errordatas)
{
    out <- dplyr::bind_rows(inside(par, errordatas[[1]]),
                            inside(par, errordatas[[2]]),
                            inside(par, errordatas[[3]]),
                            inside(par, errordatas[[4]]),
                            .id='v')
    out$v <- forcats::fct_recode(out$v,
                                 'RISK MODEL'      ='1',
                                 'TIME MODEL'      ='2',
                                 'BLOCK-DIAG MODEL'='3',
                                 'COMPLETE MODEL'  ='4')
    return(out)
}

bias2plot <- function(df, title)
{
    ggplot(df, aes(ymin=q025, ymax=q975, x=data, color=conv, size=nr))+
        geom_linerange(position=position_dodge(width=0.2))+
        geom_point(mapping=aes(x=data, y=mean), size=3,
                   shape=21, color='black', fill='white', stroke=2)+
        facet_grid(~v, scales='free_x')+
        geom_hline(yintercept=0, linetype='dashed')+
        labs(x=NULL, y='Bias',
             title=title,
             subtitle='with 2.5% and 97.5% quantiles',
             color='Convergance\n(%)',
             size='Number of\nmodels')+
        coord_flip()+
        ## theme_minimal()+
        ## scale_color_viridis_c()+
        ## scale_color_distiller(palette='Spectral')+
        scale_color_fish(option='Trimma_lantana')+
        guides(color=guide_colorbar(barwidth=20))+
        theme(
            strip.background=element_rect(colour='black', fill='white'),
            legend.position='bottom',
            plot.title=element_text(size=15), 
            plot.subtitle=element_text(size=14, face='italic'),
            axis.text.x=element_text(size=12), 
            axis.text.y=element_text(size=13),
            axis.title.x=element_text(size=13,
                                      margin=unit(c(t=3, r=0, b=0, l=0),
                                                  'mm')),
            legend.justification=c(1, 0),
            legend.title=element_text(size=13), 
            legend.text=element_text(size=12),
            strip.text.x=element_text(size=13, face='bold')
        )
}
bias <- function(par, title)
{
    out <- biasformat(
        par,
        list(v1error, v2error, v3error, v4error)
    )
    bias2plot(out, title=title)
}
bias('beta1', expression(bold('Parameter:'~beta[1])))

bias('beta2', expression(bold('Parameter:'~beta[2])))

bias('gama1', expression(bold('Parameter:'~gamma[1])))

bias('gama2', expression(bold('Parameter:'~gamma[2])))

bias('w1', expression(bold('Parameter:'~w[1])))

bias('w2', expression(bold('Parameter:'~w[2])))

par <- 'logs2_1'
out <- dplyr::bind_rows(inside(par, v1error),
                        inside(par, v3error),
                        inside(par, v4error),
                        .id='v')
out$v <- forcats::fct_recode(out$v,
                             'RISK MODEL'      ='1',
                             'BLOCK-DIAG MODEL'='2',
                             'COMPLETE MODEL'  ='3')
bias2plot(out, expression(bold('Parameter:'~log(sigma[1]^2))))

par <- 'logs2_2'
out <- dplyr::bind_rows(inside(par, v1error),
                        inside(par, v3error),
                        inside(par, v4error),
                        .id='v')
out$v <- forcats::fct_recode(out$v,
                             'RISK MODEL'      ='1',
                             'BLOCK-DIAG MODEL'='2',
                             'COMPLETE MODEL'  ='3')
bias2plot(out, expression(bold('Parameter:'~log(sigma[2]^2))))

par <- 'logs2_3'
out <- dplyr::bind_rows(inside(par, v2error),
                        inside(par, v3error),
                        inside(par, v4error),
                        .id='v')
out$v <- forcats::fct_recode(out$v,
                             'TIME MODEL'      ='1',
                             'BLOCK-DIAG MODEL'='2',
                             'COMPLETE MODEL'  ='3')
bias2plot(out, expression(bold('Parameter:'~log(sigma[3]^2))))

par <- 'logs2_4'
out <- dplyr::bind_rows(inside(par, v2error),
                        inside(par, v3error),
                        inside(par, v4error),
                        .id='v')
out$v <- forcats::fct_recode(out$v,
                             'TIME MODEL'      ='1',
                             'BLOCK-DIAG MODEL'='2',
                             'COMPLETE MODEL'  ='3')
bias2plot(out, expression(bold('Parameter:'~log(sigma[4]^2))))

par <- 'rhoZ12'
out <- dplyr::bind_rows(inside(par, v1error),
                        inside(par, v3error),
                        inside(par, v4error),
                        .id='v')
out$v <- forcats::fct_recode(out$v,
                             'RISK MODEL'      ='1',
                             'BLOCK-DIAG MODEL'='2',
                             'COMPLETE MODEL'  ='3')
bias2plot(out, expression(bold('Parameter:'~z(rho[12]))))

par <- 'rhoZ34'
out <- dplyr::bind_rows(inside(par, v2error),
                        inside(par, v3error),
                        inside(par, v4error),
                        .id='v')
out$v <- forcats::fct_recode(out$v,
                             'TIME MODEL'      ='1',
                             'BLOCK-DIAG MODEL'='2',
                             'COMPLETE MODEL'  ='3')
bias2plot(out, expression(bold('Parameter:'~z(rho[34]))))

bias2plot_rho <- function(df, title)
{
    ggplot(df, aes(ymin=q025, ymax=q975, x=data, color=conv, size=nr))+
        geom_linerange(position=position_dodge(width=0.2))+
        geom_point(mapping=aes(x=data, y=mean), size=3,
                   shape=21, color='black', fill='white', stroke=2)+
        facet_grid(~v, scales='free_x', labeller=label_parsed)+
        geom_hline(yintercept=0, linetype='dashed')+
        labs(x=NULL, y='Bias',
             title=title,
             subtitle='with 2.5% and 97.5% quantiles',
             color='Convergance\n(%)',
             size='Number of\nmodels')+
        coord_flip()+
        ## theme_minimal()+
        ## scale_color_viridis_c()+
        scale_color_fish(option='Trimma_lantana')+
        guides(color=guide_colorbar(barwidth=12.5))+
        theme(
            strip.background=element_rect(colour='black', fill='white'),
            legend.position='bottom',
            plot.title=element_text(size=15), 
            plot.subtitle=element_text(size=14, face='italic'),
            axis.text.x=element_text(size=12), 
            axis.text.y=element_text(size=13),
            axis.title.x=element_text(size=13,
                                      margin=unit(c(t=3, r=0, b=0, l=0),
                                                  'mm')),
            legend.justification=c(1, 0),
            legend.title=element_text(size=13), 
            legend.text=element_text(size=12),
            strip.text.x=element_text(size=13)
        )
}

out <- dplyr::bind_rows(inside('rhoZ13', v4error),
                        inside('rhoZ24', v4error),
                        inside('rhoZ14', v4error),
                        inside('rhoZ23', v4error),
                        .id='v')
out$v <- forcats::fct_recode(out$v,
                             'bold(z(rho[13]))'='1',
                             'bold(z(rho[24]))'='2',
                             'bold(z(rho[14]))'='3',
                             'bold(z(rho[23]))'='4')
bias2plot_rho(out,
              expression(bold("Complete model's cross-correlations")))

BIAS SD

bias2plot2 <- function(df, title)
{
    ggplot(df, aes(ymin=mean-1.96*sd, ymax=mean+1.96*sd, x=data,
                   color=conv, size=nr))+
        geom_linerange(position=position_dodge(width=0.2))+
        geom_point(mapping=aes(x=data, y=mean), size=3,
                   shape=21, color='black', fill='white', stroke=2)+
        facet_grid(~v, scales='free_x')+
        geom_hline(yintercept=0, linetype='dashed')+
        labs(x=NULL, y='Bias',
             title=title,
             subtitle='with \u00b1 1.96 standard deviations',
             color='Convergance\n(%)',
             size='Number of\nmodels')+
        coord_flip()+
        ## theme_minimal()+
        ## scale_color_viridis_c()+
        ## scale_color_distiller(palette='Greys', direction=1)+
        scale_color_fish(option='Trimma_lantana')+
        guides(color=guide_colorbar(barwidth=20))+
        theme(
            strip.background=element_rect(colour='black', fill='white'),
            legend.position='bottom',
            plot.title=element_text(size=15), 
            plot.subtitle=element_text(size=14, face='italic'),
            axis.text.x=element_text(size=12), 
            axis.text.y=element_text(size=13),
            axis.title.x=element_text(size=13,
                                      margin=unit(c(t=3, r=0, b=0, l=0),
                                                  'mm')),
            legend.justification=c(1, 0),
            legend.title=element_text(size=13), 
            legend.text=element_text(size=12),
            strip.text.x=element_text(size=13, face='bold')
        )
}
bias2 <- function(par, title)
{
    out <- biasformat(
        par,
        list(v1error, v2error, v3error, v4error)
    )
    bias2plot2(out, title=title)
}
bias2('beta1', expression(bold('Parameter:'~beta[1])))

bias2('beta2', expression(bold('Parameter:'~beta[2])))

bias2('gama1', expression(bold('Parameter:'~gamma[1])))

bias2('gama2', expression(bold('Parameter:'~gamma[2])))

bias2('w1', expression(bold('Parameter:'~w[1])))

bias2('w2', expression(bold('Parameter:'~w[2])))

par <- 'logs2_1'
out <- dplyr::bind_rows(inside(par, v1error),
                        inside(par, v3error),
                        inside(par, v4error),
                        .id='v')
out$v <- forcats::fct_recode(out$v,
                             'RISK MODEL'      ='1',
                             'BLOCK-DIAG MODEL'='2',
                             'COMPLETE MODEL'  ='3')
bias2plot2(out, expression(bold('Parameter:'~log(sigma[1]^2))))

par <- 'logs2_2'
out <- dplyr::bind_rows(inside(par, v1error),
                        inside(par, v3error),
                        inside(par, v4error),
                        .id='v')
out$v <- forcats::fct_recode(out$v,
                             'RISK MODEL'      ='1',
                             'BLOCK-DIAG MODEL'='2',
                             'COMPLETE MODEL'  ='3')
bias2plot2(out, expression(bold('Parameter:'~log(sigma[2]^2))))

par <- 'logs2_3'
out <- dplyr::bind_rows(inside(par, v2error),
                        inside(par, v3error),
                        inside(par, v4error),
                        .id='v')
out$v <- forcats::fct_recode(out$v,
                             'TIME MODEL'      ='1',
                             'BLOCK-DIAG MODEL'='2',
                             'COMPLETE MODEL'  ='3')
bias2plot2(out, expression(bold('Parameter:'~log(sigma[3]^2))))

par <- 'logs2_4'
out <- dplyr::bind_rows(inside(par, v2error),
                        inside(par, v3error),
                        inside(par, v4error),
                        .id='v')
out$v <- forcats::fct_recode(out$v,
                             'TIME MODEL'      ='1',
                             'BLOCK-DIAG MODEL'='2',
                             'COMPLETE MODEL'  ='3')
bias2plot2(out, expression(bold('Parameter:'~log(sigma[4]^2))))

par <- 'rhoZ12'
out <- dplyr::bind_rows(inside(par, v1error),
                        inside(par, v3error),
                        inside(par, v4error),
                        .id='v')
out$v <- forcats::fct_recode(out$v,
                             'RISK MODEL'      ='1',
                             'BLOCK-DIAG MODEL'='2',
                             'COMPLETE MODEL'  ='3')
bias2plot2(out, expression(bold('Parameter:'~z(rho[12]))))

par <- 'rhoZ34'
out <- dplyr::bind_rows(inside(par, v2error),
                        inside(par, v3error),
                        inside(par, v4error),
                        .id='v')
out$v <- forcats::fct_recode(out$v,
                             'TIME MODEL'      ='1',
                             'BLOCK-DIAG MODEL'='2',
                             'COMPLETE MODEL'  ='3')
bias2plot2(out, expression(bold('Parameter:'~z(rho[34]))))

bias2plot2_rho <- function(df, title)
{
    ggplot(df, aes(ymin=mean-1.96*sd, ymax=mean+1.96*sd, x=data,
                   color=conv, size=nr))+
        geom_linerange(position=position_dodge(width=0.2))+
        geom_point(mapping=aes(x=data, y=mean), size=3,
                   shape=21, color='black', fill='white', stroke=2)+
        facet_grid(~v, scales='free_x', labeller=label_parsed)+
        geom_hline(yintercept=0, linetype='dashed')+
        labs(x=NULL, y='Bias',
             title=title,
             subtitle='with \u00b1 1.96 standard deviations',
             color='Convergance\n(%)',
             size='Number of\nmodels')+
        coord_flip()+
        ## theme_minimal()+
        ## scale_color_viridis_c()+
        scale_color_fish(option='Trimma_lantana')+
        guides(color=guide_colorbar(barwidth=12.5))+
        theme(
            strip.background=element_rect(colour='black', fill='white'),
            legend.position='bottom',
            plot.title=element_text(size=15), 
            plot.subtitle=element_text(size=14, face='italic'),
            axis.text.x=element_text(size=12), 
            axis.text.y=element_text(size=13),
            axis.title.x=element_text(size=13,
                                      margin=unit(c(t=3, r=0, b=0, l=0),
                                                  'mm')),
            legend.justification=c(1, 0),
            legend.title=element_text(size=13), 
            legend.text=element_text(size=12),
            strip.text.x=element_text(size=13)
        )
}

out <- dplyr::bind_rows(inside('rhoZ13', v4error),
                        inside('rhoZ24', v4error),
                        inside('rhoZ14', v4error),
                        inside('rhoZ23', v4error),
                        .id='v')
out$v <- forcats::fct_recode(out$v,
                             'bold(z(rho[13]))'='1',
                             'bold(z(rho[24]))'='2',
                             'bold(z(rho[14]))'='3',
                             'bold(z(rho[23]))'='4')
bias2plot2_rho(out,
               expression(bold("Complete model's cross-correlations")))

CIF

cmean <- function(data, data_label)
{
    df  <- data%>%dplyr::select(beta1, beta2, gama1, gama2, w1, w2)
    out <-
        df%>%
        dplyr::summarize_all(mean)%>%
        dplyr::mutate(data=data_label)
    return(out)
}
lcmean <- function(datas, data_labels)
{
    out <- NULL
    for (i in seq(18))
    {
        out <- out%>%bind_rows(cmean(datas[[i]], data_labels[[i]]))
    }
    return(out)
}
df_cmean <-
    dplyr::bind_rows(lcmean(v1data, label), 
                     lcmean(v2data, label), 
                     lcmean(v3data, label), 
                     lcmean(v4data, label), .id='v')%>%
    dplyr::mutate(v=v%>%
                      forcats::fct_recode('RISK MODEL'      ='1',
                                          'TIME MODEL'      ='2',
                                          'BLOCK-DIAG MODEL'='3',
                                          'COMPLETE MODEL'  ='4'))
time  <- seq(30, 79.9, length.out=100)
delta <- 80

compCifs <- function(df_cmean.obj)
{
    beta1 <- unlist(df_cmean.obj%>%dplyr::select(beta1))
    n     <- length(beta1)
    beta2 <- unlist(df_cmean.obj%>%dplyr::select(beta2))
    gama1 <- unlist(df_cmean.obj%>%dplyr::select(gama1))
    gama2 <- unlist(df_cmean.obj%>%dplyr::select(gama2))
    w1    <- unlist(df_cmean.obj%>%dplyr::select(w1))
    w2    <- unlist(df_cmean.obj%>%dplyr::select(w2))
    cif1  <- vector(mode='list', length=n)
    cif2  <- vector(mode='list', length=n)
    for (i in seq(n))
    {
        risk1     <- exp(beta1[i])
        risk2     <- exp(beta2[i])
        level     <- 1+risk1+risk2
        gt        <- atanh(2*time/delta-1)
        traje1    <- pnorm(w1[i]*gt-gama1[i])
        traje2    <- pnorm(w2[i]*gt-gama2[i])
        cif1[[i]] <- risk1/level*traje1
        cif2[[i]] <- risk2/level*traje2
    }
    return(list(cif1=cif1, cif2=cif2))
}
cifsHigh <- compCifs(dplyr::filter(df_cmean, grepl('high', data)))
cifsLow  <- compCifs(dplyr::filter(df_cmean, grepl( 'low', data)))

getCIF <- function(compCifs.obj, n)
{
    cif1 <- cif2 <- NULL
    for (i in seq(n))
    {
        cif1 <- c(cif1, unlist(compCifs.obj[['cif1']][i]))
        cif2 <- c(cif2, unlist(compCifs.obj[['cif2']][i]))
    }
    return(cbind(cif1, cif2))
}
cifs_high <- getCIF(cifsHigh, 36)
cifs_low  <- getCIF(cifsLow,  36)

times <- rep(time, 36)
model <- rep(c('RISK MODEL',
               'TIME MODEL',
               'BLOCK-DIAG MODEL',
               'COMPLETE MODEL'), each=900)

labels <- unique(sub('-high|--low', '', label))
labels <- rep(rep(labels, each=100), 4)

dfcif_type <- function(data)
{
    out <- tibble::tibble(cif  =data,
                          time =times,
                          model=forcats::as_factor(model),
                          label=labels)
    return(out)
}
dfcif1_high <- dfcif_type(cifs_high[, 'cif1'])
dfcif2_high <- dfcif_type(cifs_high[, 'cif2'])

dfcif1_low  <- dfcif_type(cifs_low[, 'cif1'])
dfcif2_low  <- dfcif_type(cifs_low[, 'cif2'])

cif2plot <- function(dfcif_type, descr)
{
    ggplot(dfcif_type, aes(x=time, y=cif, color=label))+
        geom_line(size=1.5)+
        facet_grid(~model)+
        theme_minimal()+
        labs(x='Time', y=NULL, color=NULL,
             title=descr, subtitle='True curve in dashed black')+
        ## scale_color_fish_d(option='Scarus_quoyi')+
        scale_color_brewer(palette='Spectral')+
        theme(
            axis.text.x=element_text(size=13), 
            axis.text.y=element_text(size=13), 
            ## legend.position=c(0.05, 0.65),
            legend.position='top',
            legend.text=element_text(size=14),
            legend.box.background=element_rect(color='black'), 
            strip.background=element_rect(colour='black', fill='white'),
            strip.text.x=element_text(size=14, face='bold'),
            plot.title=element_text(size=15, face='bold'), 
            plot.subtitle=element_text(size=15, face='italic'),
            axis.title.x=element_text(size=14,
                                      margin=unit(c(t=3, r=0, b=0, l=0),
                                                  'mm'))
        )+
        guides(color=guide_legend(nrow=3))
}
truefixed <-
    tibble::tribble(
                ~beta1, ~beta2, ~gama1, ~gama2, ~w1, ~w2, 
                     3,    2.6,    2.5,      4,   5,  10,
                    -2,   -1.5,      1,    1.5,   3,   4, 
                )
ciftrue_high <- compCifs(truefixed[1, ])
ciftrue_low  <- compCifs(truefixed[2, ])

ciftrue <- function(data)
{
    out <- tibble::tibble(time=time, cif=unlist(data))
    return(out)
}

cif1true_high <- ciftrue(ciftrue_high[['cif1']])
cif2true_high <- ciftrue(ciftrue_high[['cif2']])

cif1true_low  <- ciftrue(ciftrue_low[['cif1']])
cif2true_low  <- ciftrue(ciftrue_low[['cif2']])

p1 <-
    cif2plot(dfcif1_high, descr='CIF of failure cause 1')+
    geom_line(cif1true_high, mapping=aes(x=time, y=cif),
              color='black', size=1.5, linetype='dashed')
p2 <-
    cif2plot(dfcif2_high, descr='CIF of failure cause 2')+
    geom_line(cif2true_high, mapping=aes(x=time, y=cif),
              color='black', size=1.5, linetype='dashed')+
    theme(legend.position='none')
p1/p2+plot_annotation(title='HIGH CIF SCENARIO',
                      theme=theme(plot.title=element_text(face='bold')))

p3 <-
    cif2plot(dfcif1_low, descr='CIF of failure cause 1')+
    geom_line(cif1true_low, mapping=aes(x=time, y=cif),
              color='black', size=1.5, linetype='dashed')
p4 <-
    cif2plot(dfcif2_low, descr='CIF of failure cause 2')+
    geom_line(cif2true_low, mapping=aes(x=time, y=cif),
              color='black', size=1.5, linetype='dashed')+
    theme(legend.position='none')
p3/p4+plot_annotation(title='LOW CIF SCENARIO',
                      theme=theme(plot.title=element_text(face='bold')))

HISTO LOGS2

data.label <- function(data, label)
{
    out <- data%>%dplyr::mutate(scenario=label)
    return(out)
}
logs2pilha <- function(data1, data2, data3, par, label)
{
    out1 <- out2 <- out3 <- NULL
    for (i in seq(18))
    {
        out1 <-
            out1%>%
            dplyr::bind_rows(data.label(data1[[i]]%>%dplyr::select(par),
                                        label[[i]]))
        out2 <-
            out2%>%
            dplyr::bind_rows(data.label(data2[[i]]%>%dplyr::select(par),
                                        label[[i]]))
        out3 <-
            out3%>%
            dplyr::bind_rows(data.label(data3[[i]]%>%dplyr::select(par),
                                        label[[i]]))
    }
    out <-
        dplyr::bind_rows(out1, out2, out3, .id='v')%>%
        dplyr::filter(scenario == 'cs02-high-60k' |
                      scenario == 'cs05-high-60k' |
                      scenario == 'cs10-high-60k')
    return(out)
}
histo <- function(data, true, title)
{
    ggplot(data, aes(x=par, fill=scenario))+
        geom_density(alpha=0.6)+
        facet_wrap(~v, scales='free')+
        scale_fill_brewer(palette='Spectral')+
        ## scale_fill_viridis_d()+
        geom_vline(xintercept=true, linetype='dashed', size=1)+
        labs(x=NULL, y=NULL, fill=NULL,
             title=title, subtitle='True value in dashed black')+
        theme_minimal()+
        theme(
            legend.position='bottom',
            legend.box.background=element_rect(color='black'),
            legend.text=element_text(size=14), 
            plot.title=element_text(size=15),
            plot.subtitle=element_text(size=13, face='italic'),
            strip.background=element_rect(color='black', fill='white'),
            strip.text.x=element_text(size=13, face='bold'),
            axis.text.x=element_text(size=12), 
            axis.text.y=element_text(size=12)
        )
}
logs2_1df <-
    logs2pilha(v1data, v3data, v4data, 'logs2_1', label)%>%
    dplyr::rename(par=logs2_1)%>%
    dplyr::mutate(v=v%>%
                      forcats::fct_recode('RISK MODEL'      ='1',
                                          'BLOCK-DIAG MODEL'='2',
                                          'COMPLETE MODEL'  ='3'))
logs2_2df <-
    logs2pilha(v1data, v3data, v4data, 'logs2_2', label)%>%
    dplyr::rename(par=logs2_2)%>%
    dplyr::mutate(v=v%>%
                      forcats::fct_recode('RISK MODEL'      ='1',
                                          'BLOCK-DIAG MODEL'='2',
                                          'COMPLETE MODEL'  ='3'))
logs2_3df <-
    logs2pilha(v2data, v3data, v4data, 'logs2_3', label)%>%
    dplyr::rename(par=logs2_3)%>%
    dplyr::mutate(v=v%>%
                      forcats::fct_recode('TIME MODEL'      ='1',
                                          'BLOCK-DIAG MODEL'='2',
                                          'COMPLETE MODEL'  ='3'))
logs2_4df <-
    logs2pilha(v2data, v3data, v4data, 'logs2_4', label)%>%
    dplyr::rename(par=logs2_4)%>%
    dplyr::mutate(v=v%>%
                      forcats::fct_recode('TIME MODEL'      ='1',
                                          'BLOCK-DIAG MODEL'='2',
                                          'COMPLETE MODEL'  ='3'))
p1 <-
    histo(logs2_1df, log(1),
          expression(bold('Parameter:'~log(sigma[1]^2))))+
    theme(legend.position='none')
p2 <-
    histo(logs2_2df, log(0.6), 
          expression(bold('Parameter:'~log(sigma[2]^2))))+
    theme(legend.position='none')
p3 <-
    histo(logs2_3df, log(0.7), 
          expression(bold('Parameter:'~log(sigma[3]^2))))+
    theme(legend.position='none')
p4 <-
    histo(logs2_4df, log(0.9), 
          expression(bold('Parameter:'~log(sigma[4]^2))))
p1/p2/p3/p4

HISTO RHOZ

crossco <- function()
{
    out <- NULL
    for (i in seq(18))
    {
        out <-
            out%>%
            dplyr::bind_rows(data.label(
                       v4data[[i]]%>%
                       dplyr::select(rhoZ13, rhoZ24, rhoZ14, rhoZ23),
                       label[[i]]
                   ))
    }
    out <- out%>%
        tidyr::pivot_longer(
                   !scenario, names_to='par', values_to='value'
               )%>%
        dplyr::mutate(par=forcats::as_factor(par)%>%
                          forcats::fct_recode(
                                       'bold(z(rho[13]))'='rhoZ13',
                                       'bold(z(rho[24]))'='rhoZ24',
                                       'bold(z(rho[14]))'='rhoZ14',
                                       'bold(z(rho[23]))'='rhoZ23'))%>%
        dplyr::filter(scenario == 'cs02-high-60k' |
                      scenario == 'cs05-high-60k' |
                      scenario == 'cs10-high-60k')
    return(out)
}
rhoZ_df  <- crossco()
rhoZtrue <-
    tibble::tibble(
                par=forcats::as_factor(levels(rhoZ_df$par)),
                value=atanh(c(-0.5, -0.4, 0.3, 0.3))
            )
rhoZ12_df <-
    logs2pilha(v1data, v3data, v4data, 'rhoZ12', label)%>%
    dplyr::rename(par=rhoZ12)%>%
    dplyr::mutate(v=v%>%
                      forcats::fct_recode('RISK MODEL'      ='1',
                                          'BLOCK-DIAG MODEL'='2',
                                          'COMPLETE MODEL'  ='3'))
rhoZ34_df <-
    logs2pilha(v2data, v3data, v4data, 'rhoZ34', label)%>%
    dplyr::rename(par=rhoZ34)%>%
    dplyr::mutate(v=v%>%
                      forcats::fct_recode('TIME MODEL'      ='1',
                                          'BLOCK-DIAG MODEL'='2',
                                          'COMPLETE MODEL'  ='3'))
p5 <-
    histo(rhoZ12_df, atanh(0.1),
          expression(bold('Parameter:'~z(rho[12]))))+
    theme(legend.position='none')
p6 <-
    histo(rhoZ34_df, atanh(0.2),
          expression(bold('Parameter:'~z(rho[34]))))+
    theme(legend.position='none')
p7 <-
    ggplot(rhoZ_df, aes(x=value, fill=scenario))+
    geom_density(alpha=0.6)+
    facet_grid(~par, scales='free', labeller=label_parsed)+
    scale_fill_brewer(palette='Spectral')+
    ## scale_fill_viridis_d()+
    geom_vline(rhoZtrue, mapping=aes(xintercept=value),
               linetype='dashed', size=1)+
    labs(x=NULL, y=NULL, fill=NULL,
         title="Complete model's cross-correlations",
         subtitle='True value in dashed black')+
    theme_minimal()+
    theme(
        legend.position='bottom',
        legend.box.background=element_rect(color='black'),
        legend.text=element_text(size=14), 
        plot.title=element_text(size=15, face='bold'),
        plot.subtitle=element_text(size=13, face='italic'),
        strip.background=element_rect(color='black', fill='white'),
        strip.text.x=element_text(size=13, face='bold'),
        axis.text.x=element_text(size=12), 
        axis.text.y=element_text(size=12)
    )
p5/p6/p7

COR

thepilha <- function()
{
    out <- NULL
    for (i in seq(18))
    {
        out <-
            out%>%
            dplyr::bind_rows(data.label(v4data[[i]]%>%
                                        dplyr::select(-c(conv, mll)),
                                        label[[i]]))
    }
    out <-
        out%>%
        dplyr::filter(scenario == 'cs10-high-60k')%>%
        dplyr::select(!scenario)
    return(out)
}
v4pars    <- thepilha()
v4cor     <- round(cor(v4pars), 1)
corlabels <- c(expression(beta[1]),
               expression(beta[2]),
               expression(gamma[1]),
               expression(gamma[2]),
               expression(w[1]),
               expression(w[2]),
               expression(log(sigma[1]^2)),
               expression(log(sigma[2]^2)),
               expression(log(sigma[3]^2)),
               expression(log(sigma[4]^2)),
               expression(z(rho[12])),
               expression(z(rho[13])),
               expression(z(rho[14])),
               expression(z(rho[23])),
               expression(z(rho[24])),
               expression(z(rho[34])))

ggcorrplot(v4cor, 
           ## type='lower', 
           ## lab=TRUE, 
           title='cs10-high-60k, complete model',
           ## p.mat=cor_pmat(v4cor), 
           ## insig='blank', 
           ggtheme=theme(
               plot.title=element_text(size=19, face='bold'), 
               axis.title.x=element_blank(),
               axis.title.y=element_blank(),
               axis.ticks=element_blank(),
               legend.title=element_blank(),
               legend.text=element_text(size=15)
           ))+
    scale_x_discrete(labels=corlabels)+
    scale_y_discrete(labels=corlabels)+
    theme(axis.text.x=element_text(size=18), 
          axis.text.y=element_text(size=18))+
    ## scale_fill_viridis_c(option='plasma', direction=-1)
    scale_fill_fish(option='Trimma_lantana')+
    guides(fill=guide_colorbar(barheight=20))